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1.  Introduction 


Atomistic  simulations,  including  molecular  dynamics  (MD),  have  become  an 
important  tool  for  researchers  in  various  disciplines  to  study  complex  physical  and 
chemical  processes  on  the  molecular  scale.1  However,  there  are  2  very  important 
factors  that  affect  the  quality  of  the  MD  results:  the  description  of  the  interatomic 
potentials  and  the  initial  configuration  of  the  system  for  the  trajectory  from  which 
relevant  thermodynamic  and  structural  averages  will  be  obtained.  The  interatomic 
potentials  are  typically  fitted  to  results  from  experiment2  or  quantum-level 
atomistic  simulations,3  such  as  density  functional  theory.4  ReaxFF5  and,  more 
recently,  ReaxFF-/g6  have  emerged  as  preferred  potentials  for  reactive  systems  and 
have  been  fitted  for  use  in  simulations  of  many  types  of  materials,  such  as 
hydrocarbons,5  energetic  materials,7  semiconductors,8  and  metal  oxides.9  A 
relevant  initial  configuration  is  usually  obtained  through  an  equilibration  trajectory 
in  which  the  system  is  relaxed  away  from  the  initial  atomic  positions  and  velocities, 
through  integration  of  the  equations  of  motion.  The  length  of  the  equilibration 
trajectory  should  be  sufficiently  long  to  converge  to  the  desired  thermodynamic 
state. 

The  original  parameterization  of  the  ReaxFF-/g  potential6  was  developed  for  use  in 
CHNO-type  energetic  materials  but  has  accurately  predicted  structural  properties 
of  the  nonenergetic  molecular  crystal  sucrose  (C12H22O11),  or  table  sugar.  Our 
simulations  of  sucrose  were  started  with  the  experimental  configuration,  but  planar 
defects  were  observed  to  nucleate  during  the  equilibration  portion  of  the  simulation. 
This  technical  note  describes  a  modification  of  the  equilibration  process  that 
produces  a  defect-free  crystal  for  sucrose  using  the  ReaxFF-Zg  potential. 

2.  Defect-Free  Equilibration 

MD  simulations  of  sucrose  were  initiated  using  the  experimental  lattice  vectors  and 
atomic  positions  at  0  GPa.10  The  space  group  of  the  sucrose  crystal  is  P2i,  with  a 
unit  cell  containing  2  molecules.11  The  simulation  cell  is  fully  periodic,  containing 
3x3x3  unit  cells,  for  a  total  of  2,430  atoms.  Figure  la  shows  the  molecular  centers 
of  mass  for  this  configuration,  projected  onto  the  vz-planc.  The  timestep  of  the 
simulations  is  0.1  fs,  with  thermostating  and  barostating  rates  of  0.1  ps  and  1.0  ps, 
respectively.  All  simulations  are  performed  using  the  LAMMPS  computer  code.12 
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Fig.  1  a)  Initial  and  b)  final  configurations  of  the  molecular  centers  of  mass  (projected  onto 
the  *z-plane)  for  Equilibration  1  using  the  experimental  unit  cell  with  the  Goddard  model  at 
OGPa 

In  a  similar  procedure  to  that  of  an  earlier  study  for  simulations  of  energetic 
molecular  crystals  using  ReaxFF-/g,13  the  equilibration  starts  with  microcanonical 
(NVE)  MD  with  velocity  rescaling  every  timestep  for  2.5  ps,  followed  by 
isothermal-isochoric  (NVT)  MD  at  300  K  for  2.5  ps  with  no  rescaling,  then 
isothermal-isobaric  (NPT)  MD  at  300  K  and  0  GPa  for  2.5  ps  during  which  the  cell 
angles  are  fixed  but  edge  lengths  are  allowed  to  change.  Finally,  isothermal- 
isostress  (NsT)  MD  is  performed  for  40  ps,  during  which  cell  angles  are  allowed  to 
change  along  with  edge  lengths.  The  first  35  ps  of  the  NsT  MD  simulation  is  still 
considered  part  of  the  equilibration;  the  last  5  ps  is  the  production  portion  of  the 
trajectories.  This  equilibration  protocol  was  sufficient  to  provide  initial  states 
converged  to  proper  thermodynamic  states  for  the  energetic  material 
cyclotrimethylene  trinitramine  (RDX). 13,14 

Structural  and  thermodynamic  averages  are  taken  from  snapshots  created  every  500 
timesteps  over  the  last  5  ps.  Figure  lb  shows  the  final  configuration  of  the 
molecular  centers  of  mass  at  the  end  of  the  trajectory  using  this  equilibration 
protocol,  denoted  “Equilibration  1”;  the  Table  below  contains  the  corresponding 
lattice  constants  and  angles.  The  lattice  constants,  a,  b,  and  c,  are  all  larger  than  the 
experimental  configuration,  resulting  in  a  density  of  1.56  g/cm3.  This  is  1.89%  less 
than  the  experimental  value  of  1.59  g/cm3.  This  indicates  that  the  initial  system  in 
Fig.  la  (i.e.,  the  experimental  configuration)  is  under  compression  when  the 
equilibration  trajectory  initiates.  During  the  NVE  portion  of  the  equilibration,  in 
which  the  volume  is  fixed  and  the  compressive  stress  cannot  be  relaxed  in  any  way 
other  than  reorientation  and/or  deformation  of  the  molecules,  this  compressive 
stress  causes  the  crystal  to  rotate  around  the  y-axis,  as  seen  from  the  dashed  blue 
reference  lines  in  Fig.  la  and  lb.  The  enforcement  of  periodic  boundaries  in  the 
simulation  results  in  a  stacking  fault,  indicated  by  the  black  dashed  line  in  Fig.  lb. 
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Table  Crystal  parameters  after  each  equilibration,  compared  to  experiment 


Experiment 

Equilibration  1 

Error3 

Equilibration  2 

Error3 

a 

10.863  A 

11.128  A 

2.44% 

11.102  A 

2.20% 

b 

8.705  A 

8.771  A 

0.76% 

8.767  A 

0.71% 

c 

7.7585  A 

7.784  A 

0.33% 

7.649  A 

-1.41% 

a 

90.0° 

89.02° 

-0.98° 

89.81° 

-0.91° 

p 

102.94° 

105.83° 

2.98° 

101.18° 

-1.77° 

y 

90.0° 

91.42° 

1.42° 

90.0° 

O 

O 

d 

3  Error  is  equal  to  the  percent  difference  for  the  distances  and  absolute  difference  for  the  angles  compared  to 
experimental  values. 


The  presence  of  the  defect  is  an  artifact  of  the  equilibration  protocol  when  using  an 
initial  highly  stressed  state.  Therefore,  it  was  necessary  to  develop  a  procedure  that 
precludes  this  artifact.  This  procedure,  denoted  “Equilibration  2”,  is  as  follows:  The 
lattice  vectors  produced  by  the  simulation  resulting  from  Equilibration  1  are  used 
in  conjunction  with  the  experimental  fractional  positions  of  each  atom  to  generate 
a  new  3x3x3  system,  as  seen  in  Fig.  2a.  This  alleviates  the  compressive  stress  on 
the  system  that  was  apparent  when  the  experimental  lattice  vectors  of  the  unit  cell 
were  used  to  construct  the  simulation  cell.  This  new  system  is  then  equilibrated  in 
the  same  manner  as  described  earlier:  2.5  ps  in  the  NVE  ensemble,  2.5  ps  in  the 
NVT  ensemble,  2.5  ps  in  the  NPT  ensemble,  and  40  ps  in  the  NsT  ensemble.  Once 
again,  structural  and  thermodynamics  averages  are  obtained  every  500  timesteps 
over  the  final  5  ps  of  the  NsT  portion  of  the  simulation.  We  stress  that  this  trajectory 
is  not  a  continuation  of  the  first  equilibration  but  an  entirely  new  simulation  using 
an  initial  configuration  obtained  from  Equilibration  1.  The  molecular  centers  of 
mass  at  the  end  of  this  trajectory  are  shown  in  Fig.  2b,  and  we  note  the  absence  of 
any  defects  in  the  system.  The  lattice  constants  and  angles  appear  in  the  Table  under 
the  column  “Equilibration  2”.  The  errors  are  reduced  for  all  crystallographic 
parameters  except  the  c  lattice  constant.  Surprisingly,  the  density  after  the  second 
equilibration  is  1.56  g/cm3,  identical  to  that  of  the  first  equilibration,  but  this  is  due 
to  offsets  in  lattice  constants  and  cell  angles. 


Fig.  2  Same  as  Fig.  1,  with  a)  the  initial  configuration  generated  by  combining  the 
experimental  atomic  fractional  coordinates  and  the  lattice  vectors  from  the  first  equilibration, 
resulting  in  b)  a  defect-free  final  configuration  after  Equilibration  2 
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To  ensure  that  the  defect  resulting  from  Equilibration  1  is  not  an  artifact  of  the  small 
system  size,  the  calculation  using  Equilibration  1  was  repeated  onal5xl5xl5 
unit  cell  system  (303,750  atoms).  The  initial  molecular  centers  of  mass  are  shown 
in  Fig.  3a,  while  the  final  configuration  after  Equilibration  1  appears  in  Fig.  3b.  It 
is  clear  that  defects  were  generated  during  the  equilibration.  In  order  to  perform 
Equilibration  2,  the  average  unit  cell  lattice  vectors  were  obtained  from  the  final 
configuration  of  the  trajectory  using  Equilibration  1  (Fig.  3b).  These  were  used, 
along  with  the  experimental  fractional  atomic  positions,  to  generate  a  new  15  x  15 
x  15  unit  cell  system.  The  equilibration  trajectory  was  then  repeated.  The  result  of 
the  trajectory  using  Equilibration  2  is  in  Fig.  3c.  Just  as  in  the  smaller  system, 
Equilibration  1  results  in  a  defect  system,  while  Equilibration  2  is  a  perfect  single 
crystal.  This  suggests  that  the  phenomenon  was  not  due  to  size  effects  and  that  a 
double  equilibration  similar  to  the  process  described  herein  should  be  performed 
for  simulations  of  any  size. 


t|»**,«*#«,*#»,»*S*»,***%»* 


(b) 


Fig.  3  Same  as  previous  figures,  for  a  larger  system  (15  x  15  x  15  unit  cells):  a)  initial 
configuration,  b)  configuration  after  Equilibration  1,  and  c)  final  configuration  after 
Equilibration  2 

This  process  points  to  the  precarious  nature  of  the  ReaxFF  potential  (which  includes 
the  ReaxFF-/g  formulation)  wherein  great  care  must  be  taken  during  equilibration 
to  avoid  undesirable  effects  in  the  system  or  from  the  subsequent  portion  of  the 
simulation  in  which  averages  are  taken.  This  especially  holds  true  for  studies  such 
as  this,  where  the  potential  was  not  fitted  for  the  material  of  interest — in  this  case, 
sucrose.  In  another  study,  we  have  observed  the  nucleation  of  a  defect  region  in  the 
energetic  material  1,1 -diamino-2, 2-dinitroethene  (FOX-7)  using  a  different 
parameterization  of  ReaxFF-/g  when  fitted  to  RDX  information.13  Once  again,  we 
are  able  to  eliminate  the  FOX-7  defect  using  the  same  double  equilibration  method 
described  here.  The  nucleation  of  defects  during  equilibration  does  not  necessarily 
bring  into  question  the  transferability  of  ReaxFF  or  ReaxFF-/g,  only  that  care  must 
be  taken  in  setting  up  the  simulation  cell  used  to  initiate  the  trajectories  from  which 
relevant  structural  and  thermodynamic  properties  are  obtained. 
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3.  Summary  and  Conclusion 


We  have  developed  an  equilibration  protocol  for  ReaxFF-type  potentials  that  is 
capable  of  avoiding  the  nucleation  of  defects  due  to  highly  stressed  initial  states. 
The  equilibration  technique  is  as  follows:  A  first  equilibration  (that  results  in  a 
system  with  one  or  more  defects)  provides  lattice  vectors  that  can  be  used  to  build 
a  new  simulation  cell.  The  trajectory  is  repeated  using  the  new  simulation  cell  as 
the  initial  configuration,  resulting  in  a  defect-free  system.  This  technique  has  been 
successfully  applied  to  sucrose  and  FOX-7.  We  have  also  shown  that  the  nucleation 
of  defects  occurs  for  larger  systems,  reducing  the  likelihood  that  a  size  effect  is 
responsible  for  this  phenomenon. 

There  are  numerous  parameterizations  of  ReaxFF  and  ReaxFF-/g,  and  this 
technique  is  likely  not  necessary  for  simulations  of  systems  whose  thermodynamic 
and  structural  states  are  similar  to  information  used  in  the  fitting.  In  those  cases,  the 
parameterization  is  expected  to  result  in  simulation  cells  that  are  not  highly  stressed. 
But  the  ReaxFF  force  fields  used  in  this  work  were  not  parameterized  for  either 
sucrose  or  FOX-7,  and  the  compressive  stresses  at  the  outset  of  the  simulations 
result  in  the  nucleation  of  planar  defects.  This  technical  note  clearly  shows  that 
some  care  must  be  taken  during  equilibration  to  get  a  suitable  starting  configuration 
and  provides  a  straightforward  means  of  relaxing  systems  initially  under  stress. 
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RDRL  WML 
M  ZOLTOSKI 
RDRL  WML  B 
I BATYREV 
J  BRENNAN 
E  BYRD 
S  IZVYEKOV 
W  MATTSON 
N  S  WEINGARTEN 
BRICE 
D  TAYLOR 
N  TRIVEDI 
RDRL  WMM  F 
M  TSCHOPP 
RDRL  WMM  G 
J  ANDZELM 
T  CHANTAWANSRI 
B  RINDERSPACHER 
TSIRK 


1  CEA 

(PDF)  J-B  MAILLET 

2  DEFENCE  RSRCH  AND  DEV 
(PDF)  CANADA 

H  ABOU-RACHID 
A  HU 


1  CFD  RSRCH  CORP 
(PDF)  D  SENGUPTA 

2  RENSSELAER  POLYTECHNIC 
(PDF)  INSTITUTE 

SDE 

CPICU 


1  UNIV  OF  ILLINOIS  AT 

(PDF)  URBANA-CHAMPAIGN 
S  CAHUDHURI 


1  UNIV  OF  MARYLAND 

(PDF)  P  CHUNG 


8 


